## Visualize experiment from hives 7-8
## Bees were given two treatments -- reward first or second
#install packages
ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", 
              "plyr", "effects", "viridis", "GGally", "survminer", 
              "tidyr","gamm4")
ipak(packages)
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: sjPlot
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: plyr
## Loading required package: effects
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Loading required package: viridis
## Loading required package: viridisLite
## Loading required package: GGally
## Loading required package: survminer
## Loading required package: ggpubr
## Loading required package: magrittr
## Loading required package: tidyr
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
## The following object is masked from 'package:Matrix':
## 
##     expand
## The following object is masked from 'package:reshape2':
## 
##     smiths
## Loading required package: gamm4
## Loading required package: mgcv
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
## This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.
## This is gamm4 0.2-5
##   ggplot2  reshape2      lme4    sjPlot  multcomp      plyr   effects 
##      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE 
##   viridis    GGally survminer     tidyr     gamm4 
##      TRUE      TRUE      TRUE      TRUE      TRUE
# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))

# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"


# print system info
print(paste("last run ", Sys.time()))
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'
## [1] "last run  2017-12-09 06:53:18"
print(R.version)
##                _                           
## platform       x86_64-apple-darwin15.6.0   
## arch           x86_64                      
## os             darwin15.6.0                
## system         x86_64, darwin15.6.0        
## status                                     
## major          3                           
## minor          4.2                         
## year           2017                        
## month          09                          
## day            28                          
## svn rev        73368                       
## language       R                           
## version.string R version 3.4.2 (2017-09-28)
## nickname       Short Summer

Read in data

bees <- read.csv(file.path(dataDir, "03_StartStop_cleaned.csv"))
bees$hive = as.factor(bees$hive)
head(bees)
##   freq     amp                  datetime visitNum treatment_rewarded
## 1  370 0.79858  2017_05_30__09_52_48_207        1                  F
## 2  340 0.68663  2017_05_30__09_52_50_473        2                  F
## 3  330 0.34916  2017_05_30__09_52_52_817        3                  F
## 4  320 0.66814  2017_05_30__09_52_56_614        4                  F
## 5  310 0.90015  2017_05_30__09_52_58_795        5                  F
## 6  320 1.10685  2017_05_30__09_52_59_487        6                  F
##                     BeeNumCol                                  accFile
## 1 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_48_207_220_450_0_5.txt
## 2 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_50_473_220_450_0_5.txt
## 3 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_52_817_220_450_0_5.txt
## 4 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_56_614_220_450_0_5.txt
## 5 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_58_795_220_450_0_5.txt
## 6 B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_59_487_220_450_0_5.txt
##   beeID hive treatment totalRewards meanFreqRewarded
## 1   207    7    RewSec           50            303.8
## 2   207    7    RewSec           50            303.8
## 3   207    7    RewSec           50            303.8
## 4   207    7    RewSec           50            303.8
## 5   207    7    RewSec           50            303.8
## 6   207    7    RewSec           50            303.8
##   totalBuzzesNotRewarded meanFreqUnrewarded IT_mm   amp_acc
## 1                     50              307.2  4.11  78.52311
## 2                     50              307.2  4.11  67.51524
## 3                     50              307.2  4.11  34.33235
## 4                     50              307.2  4.11  65.69715
## 5                     50              307.2  4.11  88.51032
## 6                     50              307.2  4.11 108.83481
# calculate mean freq rewarded vs. unrewarded for each bee
# to make sure it agrees with what is already in the dataframe
tapply(bees$freq, INDEX = list(bees$beeID, bees$treatment_rewarded), mean)
##            F        T
## 207 307.2000 303.8000
## 208 344.6000 338.0000
## 209 317.8000 334.6000
## 210 342.5000 343.6000
## 211 376.3265       NA
## 212 331.5152       NA
## 213       NA 337.2093
## 214 327.8000 349.8000
## 215       NA 356.6667
## 216 366.8000 342.6000
## 217 310.0000       NA
## 218 363.2000 327.0000
## 219 353.6000 312.6000
## 220 340.4000 336.4000
## 221 358.1818       NA
## 222 331.0000 304.4000
## 223 347.5000       NA
## 224 350.8000 348.8000
## 225 347.7419 354.4000
## 226       NA 361.4634
## 227 307.4000 309.3878
## 228 336.2000 361.4000
## 229 361.2000 351.0000
## 230 336.2000 298.2000
## 231 347.8000 367.4000
## 232 358.0000 343.8000
## 233 363.5714       NA
## 234 280.0000 313.0000
## 235       NA 352.0588
## 236 394.4828       NA
## 237 333.2000 313.2000
## 238 305.3333       NA
## 239 401.4000 395.0000
## 240 342.3810       NA
## 241 343.2000 319.8000
## 242 303.1579       NA
## 243 346.2000 348.0000
## 244 380.0000 352.2000
## 245 296.3636       NA
## 246 352.4000 329.1429
## 247 314.0000 304.6000
## 248 389.2000 349.8000
## 249 346.8000 336.4000
## 250 355.6000 347.4000
## 251 364.6341       NA
## 252 299.3333       NA
## 253 380.0000       NA
## 254 321.3636 314.0000
## 255 347.2000       NA
## 256 273.7500       NA
## 257 320.0000 347.8000
## 258 302.2222       NA
## 259 330.0000       NA
## 260       NA 328.2143
## 261 328.2353       NA
## 262 343.3333       NA
## 263 336.8000 315.2000
## 264       NA 294.0000
## 265       NA 311.3043
## 266       NA 356.6667
## 267 334.5714 330.4000
## 268 328.9362       NA
## 269 329.6000 363.6000
## 270 339.6000 344.2000
## 271 333.2609       NA
## 272 327.6471       NA
## 273 357.6471 357.4000
## 274 275.5556       NA
## 275 345.6000 336.6000
## 276       NA 336.4286
## 277 353.8000 333.2000
## 278 319.0000 305.2000
## 279 307.6000 287.8000
## 280 295.0000       NA
## 281       NA 344.2105
## 282 355.2000 360.0000
## 283 356.1905       NA
## 284 323.8462       NA
## 285 345.2000 340.8000
## 286 305.4000 297.2000
## 287       NA 386.2857
## 288 348.0000 345.1724
## 289 340.9091       NA
## 290 386.4000 359.0000
## 291 401.4706 369.8000
## 292       NA 347.7273
## 293 368.9474       NA
## 294 368.6000 365.8000
## 295 331.0526       NA
## 296 316.6000 312.4000
## 297 411.7778       NA
## 298 328.4000 314.2000
## 299       NA 318.8889
## 300 297.4000 289.8000
## 301       NA 385.8824
## 302 381.2903       NA
# number of bees -- should be 96
length(unique(bees$BeeNumCol))
## [1] 96
# num observations -- should be 5851
nrow(bees)
## [1] 5851

Visualize data

ggplot(bees, aes(y = freq, x = treatment_rewarded)) +
     geom_boxplot() +  
     stat_smooth(method = "loess", span = 0.75)

bees$trt2 <- mapvalues(bees$treatment, from = c("RewFir", "RewSec"), 
                       to = c("Rewarded 1st 50 Buzzes", "Rewarded 2nd 50 Buzzes"))

ggplot(bees, aes(x = visitNum, y = freq, color = treatment_rewarded)) +
     stat_smooth(method = "loess", span = 1) + 
     xlab("Buzz Number") + 
     ylab("Buzz Frequency(hz)") + 
     facet_wrap(~trt2) + 
     scale_color_viridis(name = "Rewarded Buzzes", 
                         discrete = TRUE, begin = 0.2, end = 0.6)

#ggsave("beeRewards.png", width = 8, height = 5)


# make survival plot
survBees <- t(sapply(unique(bees$beeID), function(x) {
     tmp <- bees[bees$beeID == x, ]
     return(tmp[which.max(tmp$visitNum), ])
}))


lastObs <- as.data.frame(survBees)
lastObs$index <- as.numeric(lastObs$visitNum)


# 1 = censored, 2 = stopped buzzing for 5 min
lastObs$status <- mapvalues(lastObs$index == 100, from = c(TRUE, FALSE), to = c(1, 2))

# create survival object
lastObs$survObj <- with(lastObs, Surv(index, status == 2))
head(lastObs)
##   freq     amp datetime visitNum treatment_rewarded BeeNumCol accFile
## 1  300 0.88263      100      100                  2         1     100
## 2  320 0.90242      200      100                  2         2     200
## 3  350 0.28458      300      100                  1         3     300
## 4  310 0.84888      374       74                  1         4     374
## 5  400 0.46176      423       49                  1         5     423
## 6  330 0.57663      456       33                  1         6     456
##   beeID hive treatment totalRewards meanFreqRewarded
## 1   207    1         2           50            303.8
## 2   208    1         2           50              338
## 3   209    1         1           50            334.6
## 4   210    1         1           50            343.6
## 5   211    2         2            0               NA
## 6   212    2         2            0               NA
##   totalBuzzesNotRewarded meanFreqUnrewarded IT_mm  amp_acc trt2 index
## 1                     50              307.2  4.11 86.78761    2   100
## 2                     50              344.6  4.46 88.73353    2   100
## 3                     50              317.8   4.4  27.9823    1   100
## 4                     24              342.5  4.28 83.46903    1    74
## 5                     49           376.3265  4.96 45.40413    2    49
## 6                     33           331.5152  4.39 56.69912    2    33
##   status survObj
## 1      1    100+
## 2      1    100+
## 3      1    100+
## 4      2     74 
## 5      2     49 
## 6      2     33
lastObs$treatment <- as.factor(unlist(lastObs$treatment))


## Kaplan-Meier estimator. 
km.as.one <- survfit(survObj ~ 1, data = lastObs, conf.type = "log-log")
km.by.trt <- survfit(survObj ~ treatment, data = lastObs, conf.type = "log-log")

par(mfrow = c(1,1))
plot(km.by.trt)

gsv <- ggsurvplot(km.by.trt, CI = FALSE) 

gsv$plot + facet_wrap(~treatment) + 
     theme_bw()   + 
     theme(legend.position = "none")  + 
     scale_color_viridis(discrete = TRUE, begin = 0.3, end = 0.7, option = "magma")  + 
     xlab("Buzz Number") + 
     ylab("Proportion Sonicating") 

#ggsave("beeSonicationPercentage.pdf", width = 8, height = 5)


ggsurv(km.by.trt, CI = FALSE, surv.col = c("#a6cee3", "#1f78b4")) + 
     ylim(c(0, 1)) + 
     xlab("Buzz Number") + 
     ylab("Proportion Sonicating") 
## Loading required package: scales
## 
## Attaching package: 'scales'
## The following object is masked from 'package:viridis':
## 
##     viridis_pal

xtabs(~lastObs$treatment + as.character(lastObs$hive))
##                  as.character(lastObs$hive)
## lastObs$treatment  1  2
##            RewFir 30 20
##            RewSec 23 23

Visualize data – averaged by bee

This is not super informative

beeAvgs = bees[bees$visitNum == 1, ]
nrow(beeAvgs) # should be 96
## [1] 96
head(beeAvgs)
##     freq     amp                  datetime visitNum treatment_rewarded
## 1    370 0.79858  2017_05_30__09_52_48_207        1                  F
## 101  350 1.25158  2017_05_30__10_15_55_623        1                  F
## 201  390 1.19220  2017_05_30__10_28_20_624        1                  T
## 301  350 0.79177  2017_05_30__10_38_37_411        1                  T
## 375  390 0.78658  2017_05_30__11_13_02_616        1                  F
## 424  270 0.04932  2017_05_30__11_27_40_579        1                  F
##                       BeeNumCol                                  accFile
## 1   B_207_2017_05_30_H_7_RewSec 2017_05_30__09_52_48_207_220_450_0_5.txt
## 101 B_208_2017_05_30_H_7_RewSec 2017_05_30__10_15_55_623_220_450_0_5.txt
## 201 B_209_2017_05_30_H_7_RewFir 2017_05_30__10_28_20_624_220_450_0_5.txt
## 301 B_210_2017_05_30_H_7_RewFir 2017_05_30__10_38_37_411_220_450_0_5.txt
## 375 B_211_2017_05_30_H_8_RewSec 2017_05_30__11_13_02_616_220_450_0_5.txt
## 424 B_212_2017_05_30_H_8_RewSec 2017_05_30__11_27_40_579_220_450_0_5.txt
##     beeID hive treatment totalRewards meanFreqRewarded
## 1     207    7    RewSec           50            303.8
## 101   208    7    RewSec           50            338.0
## 201   209    7    RewFir           50            334.6
## 301   210    7    RewFir           50            343.6
## 375   211    8    RewSec            0               NA
## 424   212    8    RewSec            0               NA
##     totalBuzzesNotRewarded meanFreqUnrewarded IT_mm    amp_acc
## 1                       50           307.2000  4.11  78.523107
## 101                     50           344.6000  4.46 123.065880
## 201                     50           317.8000  4.40 117.227139
## 301                     24           342.5000  4.28  77.853491
## 375                     49           376.3265  4.96  77.343166
## 424                     33           331.5152  4.39   4.849558
##                       trt2
## 1   Rewarded 2nd 50 Buzzes
## 101 Rewarded 2nd 50 Buzzes
## 201 Rewarded 1st 50 Buzzes
## 301 Rewarded 1st 50 Buzzes
## 375 Rewarded 2nd 50 Buzzes
## 424 Rewarded 2nd 50 Buzzes
data_long <- gather(beeAvgs[, c("BeeNumCol", "meanFreqRewarded", "meanFreqUnrewarded")], condition, meanFreq, meanFreqRewarded, meanFreqUnrewarded)
head(data_long)
##                     BeeNumCol        condition meanFreq
## 1 B_207_2017_05_30_H_7_RewSec meanFreqRewarded    303.8
## 2 B_208_2017_05_30_H_7_RewSec meanFreqRewarded    338.0
## 3 B_209_2017_05_30_H_7_RewFir meanFreqRewarded    334.6
## 4 B_210_2017_05_30_H_7_RewFir meanFreqRewarded    343.6
## 5 B_211_2017_05_30_H_8_RewSec meanFreqRewarded       NA
## 6 B_212_2017_05_30_H_8_RewSec meanFreqRewarded       NA
ggplot(data_long, aes(x = condition, y = meanFreq)) + 
  geom_boxplot()
## Warning: Removed 46 rows containing non-finite values (stat_boxplot).

# look at diffs
beeAvgs$buzzDiff = beeAvgs$meanFreqRewarded - beeAvgs$meanFreqUnrewarded
hist(beeAvgs$buzzDiff)

Analysis with start/stop trials

GAMM4 shows time series

# start with gamm so I can show change by visit number
g01 = gamm4(freq ~ s(visitNum, by = trt2) + IT_mm + hive , random =  ~(1|beeID), data = bees)
par(mfrow = c(2,2))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)

summary(g01$gam) # Summary for paper 
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## freq ~ s(visitNum, by = trt2) + IT_mm + hive
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  401.540     27.410  14.649  < 2e-16 ***
## IT_mm        -16.875      6.467  -2.609   0.0091 ** 
## hive8         22.122      5.087   4.349 1.39e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                          edf Ref.df      F p-value    
## s(visitNum):trt2Rewarded 1st 50 Buzzes 4.550  4.550 28.068  <2e-16 ***
## s(visitNum):trt2Rewarded 2nd 50 Buzzes 1.332  1.332  1.267   0.359    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0627   
## lmer.REML =  59092  Scale est. = 1358.7    n = 5851
summary(g01$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 59092
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8646 -0.4554  0.1523  0.6384  3.5375 
## 
## Random effects:
##  Groups   Name                                   Variance Std.Dev.
##  beeID    (Intercept)                             551.509 23.484  
##  Xr.0     s(visitNum):trt2Rewarded 2nd 50 Buzzes    7.599  2.757  
##  Xr       s(visitNum):trt2Rewarded 1st 50 Buzzes  657.613 25.644  
##  Residual                                        1358.685 36.860  
## Number of obs: 5851, groups:  beeID, 96; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## X(Intercept)                               401.5401    27.4100  14.649
## XIT_mm                                     -16.8750     6.4675  -2.609
## Xhive8                                      22.1223     5.0873   4.349
## Xs(visitNum):trt2Rewarded 1st 50 BuzzesFx1  -7.6413     5.6183  -1.360
## Xs(visitNum):trt2Rewarded 2nd 50 BuzzesFx1   0.8391     1.4259   0.588
## 
## Correlation of Fixed Effects:
##             X(Int) XIT_mm Xhive8 X(N15B
## XIT_mm      -0.993                     
## Xhive8       0.097 -0.176              
## X(N):2R150B  0.004 -0.004  0.006       
## X(N):2R250B  0.002  0.001  0.037  0.001
dev.off()
## null device 
##           1
# save gamm plots
pdf(file.path(figDir, "Gamm_startStop_1st50.pdf"), width = 4, height = 3)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " frequency (Hz)"))
plot.gam(select = 1, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded first", side = 2, padj = -3.5)
abline(v=50, lty = 2)
dev.off()
## null device 
##           1
pdf(file.path(figDir, "Gamm_startStop_2nd50.pdf"), width = 4, height = 3)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " frequency (Hz)"))
plot.gam(select = 2, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded second", side = 2, padj = -3.5)
dev.off()
## null device 
##           1

LMER analysis

refref: todo – more model selection and diagnostics

– get p-vals, maybe?

mod1 <- lmer(freq ~ trt2 * treatment_rewarded + IT_mm + hive + (1|beeID), data = bees)
summary(mod1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 * treatment_rewarded + IT_mm + hive + (1 | beeID)
##    Data: bees
## 
## REML criterion at convergence: 59109.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8517 -0.4529  0.1408  0.6340  3.5575 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  540.1   23.24   
##  Residual             1367.5   36.98   
## Number of obs: 5851, groups:  beeID, 96
## 
## Fixed effects:
##                                                 Estimate Std. Error
## (Intercept)                                      418.664     27.882
## trt2Rewarded 2nd 50 Buzzes                       -16.160      5.091
## treatment_rewarded T                             -12.449      1.301
## IT_mm                                            -18.335      6.479
## hive8                                             23.248      5.074
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T   14.994      2.404
##                                                 t value
## (Intercept)                                      15.016
## trt2Rewarded 2nd 50 Buzzes                       -3.174
## treatment_rewarded T                             -9.566
## IT_mm                                            -2.830
## hive8                                             4.582
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T   6.238
## 
## Correlation of Fixed Effects:
##             (Intr) t2R25B trtm_T IT_mm  hive8 
## trt2Rwr250B -0.227                            
## trtmnt_rwrT -0.046  0.176                     
## IT_mm       -0.990  0.150  0.017              
## hive8        0.122 -0.131 -0.024 -0.192       
## tr2R250B:_T  0.012 -0.147 -0.542  0.002  0.041
plot(mod1)

# diagnostics -- use REML = TRUE
m1 <- update(mod1, .~., REML =TRUE)
summary(m1) # summary for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | beeID) +  
##     trt2:treatment_rewarded
##    Data: bees
## 
## REML criterion at convergence: 59109.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8517 -0.4529  0.1408  0.6340  3.5575 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  540.1   23.24   
##  Residual             1367.5   36.98   
## Number of obs: 5851, groups:  beeID, 96
## 
## Fixed effects:
##                                                 Estimate Std. Error
## (Intercept)                                      418.664     27.882
## trt2Rewarded 2nd 50 Buzzes                       -16.160      5.091
## treatment_rewarded T                             -12.449      1.301
## IT_mm                                            -18.335      6.479
## hive8                                             23.248      5.074
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T   14.994      2.404
##                                                 t value
## (Intercept)                                      15.016
## trt2Rewarded 2nd 50 Buzzes                       -3.174
## treatment_rewarded T                             -9.566
## IT_mm                                            -2.830
## hive8                                             4.582
## trt2Rewarded 2nd 50 Buzzes:treatment_rewarded T   6.238
## 
## Correlation of Fixed Effects:
##             (Intr) t2R25B trtm_T IT_mm  hive8 
## trt2Rwr250B -0.227                            
## trtmnt_rwrT -0.046  0.176                     
## IT_mm       -0.990  0.150  0.017              
## hive8        0.122 -0.131 -0.024 -0.192       
## tr2R250B:_T  0.012 -0.147 -0.542  0.002  0.041
# pval for paper
m22 <- update(m1, .~. - trt2 : treatment_rewarded , REML = TRUE)
anova(m1, m22) # p-val for interaction
## refitting model(s) with ML (instead of REML)
## Data: bees
## Models:
## m22: freq ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | beeID)
## m1: freq ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | beeID) + 
## m1:     trt2:treatment_rewarded
##     Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m22  7 59187 59234 -29587    59173                             
## m1   8 59150 59204 -29567    59134 38.864      1  4.544e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(m1)

qqnorm(ranef(m1)$beeID[[1]])
qqline(ranef(m1)$beeID[[1]])

sjp.lmer(m1, type = "re", sort = TRUE) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(m1, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## `sjp.lmer()` will become deprecated in the future. Please use `plot_model()` instead.
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

Bootstrap CI’s for figure for paper

# set number of bootstrap replicates for models
nbootSims = 10

table(bees$hive) # more trials from hive 3
## 
##    7    8 
## 3887 1964
# using hive 7, since it's the one with the most data

# calculate an average IT for prediction
ITmean = mean(tapply(bees$IT_mm, INDEX = bees$beeID, FUN = function(x) x[1] ))

pframe <- data.frame(trt2 = rep(levels(droplevels(bees$trt2)), 2),
                     treatment_rewarded = rep(levels(bees$treatment_rewarded), each = 2),
                     IT_mm = ITmean, 
                     hive = factor(7, levels = levels(bees$hive)),  
                     beeID = 99999)
pframe$freq <- 0
pp <- predict(m1, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0



### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m1, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 10"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe
##                     trt2 treatment_rewarded    IT_mm hive beeID freq
## 1 Rewarded 1st 50 Buzzes                  F 4.267292    7 99999    0
## 2 Rewarded 2nd 50 Buzzes                  F 4.267292    7 99999    0
## 3 Rewarded 1st 50 Buzzes                  T 4.267292    7 99999    0
## 4 Rewarded 2nd 50 Buzzes                  T 4.267292    7 99999    0
##        blo      bhi predMean
## 1 332.3320 345.6363 340.4257
## 2 315.8471 328.5629 324.2655
## 3 322.0536 333.3766 327.9763
## 4 318.6930 331.4843 326.8099

Make frequency plots for paper

# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
facLevs = c("Rewarded 1st 50 Buzzes\n (buzz #51-100)", 
               "Rewarded 2nd 50 Buzzes\n (buzz #1-50)", 
               "Rewarded 1st 50 Buzzes\n (buzz #1-50)",
               "Rewarded 2nd 50 Buzzes\n (buzz #50-100)"
               )

pframe$t2 <- factor(facLevs,  levels = facLevs[c(3,1,2,4)])


g00 <- ggplot(pframe, aes(x=t2, color = treatment_rewarded, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Treatment") + 
    theme(legend.position = c(.74,.8), 
          legend.background = element_rect(fill=alpha('gray90', 0.5)), 
          axis.text.x = element_text(angle = 25, hjust = 0.9)) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded?", begin =0.3, end = 0.8) + 
  geom_vline(aes(xintercept = 2.5))
g00

ggsave(plot = g00, filename = file.path(figDir, "StartStopPreds_freq2.pdf"), width =5, height = 3.5)


fac2Levs = c("51-100","1-50", "1-50","51-100")
              

pframe$buzzNums <- factor(fac2Levs,  levels = fac2Levs[c(3,1)])
pframe$beesRewarded = mapvalues(pframe$treatment_rewarded, from  = c(" F", " T"), 
                                to = c("No", "Yes"))

g001 <- ggplot(pframe, aes(x=buzzNums, color = beesRewarded, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Sonication Number") + 
    theme(legend.position = c(.74,.8), 
          legend.background = element_rect(fill=alpha('gray90', 0.5))) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded", begin =0.3, end = 0.8) +
  facet_wrap(~trt2) + 
  ylim(c(315, 357))
g001

ggsave(plot = g001, filename = file.path(figDir, "StartStopPreds_freq3.pdf"), width =5, height = 3.5)

g0 <- ggplot(pframe, aes(x=trt2, color = treatment_rewarded, y=predMean))+
     geom_point(position = position_dodge(width = 0.4))+ 
     labs(y = "Sonication frequency (Hz)", x = "Treatment") + 
    theme(legend.position = c(.74,.8), 
          legend.background = element_rect(fill=alpha('gray90', 0.5))) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1, position = position_dodge(width = 0.4)) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded", begin =0.3, end = 0.8) + 
  geom_vline(aes(xintercept = 1.5))
g0

ggsave(plot = g0, filename = file.path(figDir, "StartStopPreds_freq.pdf"), width =4, height = 3)

Amplitude analysis for paper

mod0 <- lmer(log(amp_acc) ~ trt2 * treatment_rewarded + IT_mm + hive + (1|beeID), data = bees, REML = FALSE)
mod1 <- lmer(log(amp_acc) ~ trt2+ treatment_rewarded + IT_mm + hive + (1|beeID), data = bees, REML = FALSE)

BIC(mod1, mod0) # keep mod1
##      df      BIC
## mod1  7 11347.39
## mod0  8 11354.21
anova(mod0, mod1) #keep mod1
## Data: bees
## Models:
## mod1: log(amp_acc) ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | 
## mod1:     beeID)
## mod0: log(amp_acc) ~ trt2 * treatment_rewarded + IT_mm + hive + (1 | 
## mod0:     beeID)
##      Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod1  7 11301 11347 -5643.3    11287                         
## mod0  8 11301 11354 -5642.4    11285 1.8507      1     0.1737
mod2 <- update(mod1, .~. - trt2)
BIC(mod1, mod2) # keep mod2
##      df      BIC
## mod1  7 11347.39
## mod2  6 11342.66
anova(mod1, mod2) # disagrees with BIC
## Data: bees
## Models:
## mod2: log(amp_acc) ~ treatment_rewarded + IT_mm + hive + (1 | beeID)
## mod1: log(amp_acc) ~ trt2 + treatment_rewarded + IT_mm + hive + (1 | 
## mod1:     beeID)
##      Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mod2  6 11303 11343 -5645.3    11291                           
## mod1  7 11301 11347 -5643.3    11287 3.9494      1    0.04689 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- update(mod2, .~., REML = TRUE)


summary(m2) # summary for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ treatment_rewarded + IT_mm + hive + (1 | beeID)
##    Data: bees
## 
## REML criterion at convergence: 11310.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5886 -0.4828  0.1682  0.6921  2.4185 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.04746  0.2179  
##  Residual             0.39082  0.6252  
## Number of obs: 5851, groups:  beeID, 96
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           2.41964    0.26911   8.991
## treatment_rewarded T -0.12117    0.01804  -6.718
## IT_mm                 0.35451    0.06346   5.586
## hive8                -0.25227    0.05043  -5.002
## 
## Correlation of Fixed Effects:
##             (Intr) trtm_T IT_mm 
## trtmnt_rwrT -0.035              
## IT_mm       -0.992  0.004       
## hive8        0.103  0.024 -0.182
# pval for paper
m22 <- update(m2, .~. - treatment_rewarded , REML = TRUE)
anova(m2, m22) # p-val for interaction
## refitting model(s) with ML (instead of REML)
## Data: bees
## Models:
## m22: log(amp_acc) ~ IT_mm + hive + (1 | beeID)
## m2: log(amp_acc) ~ treatment_rewarded + IT_mm + hive + (1 | beeID)
##     Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m22  5 11346 11379 -5667.8    11336                             
## m2   6 11303 11343 -5645.3    11291 44.929      1  2.043e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(m2)

qqnorm(ranef(m2)$beeID[[1]])
qqline(ranef(m2)$beeID[[1]])

sjp.lmer(m2, type = "re", sort = TRUE) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(m2, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

Bootstrap Amp CI’s for figure for paper

# set number of bootstrap replicates for models
nbootSims = 100

table(bees$hive) # more trials from hive 3
## 
##    7    8 
## 3887 1964
# using hive 7, since it's the one with the most data

pframe <- data.frame(trt2 = rep(levels(droplevels(bees$trt2)), 2),
                     treatment_rewarded = rep(levels(bees$treatment_rewarded), each = 2),
                     IT_mm = ITmean, 
                     hive = factor(7, levels = levels(bees$hive)),  
                     beeID = 99999)
pframe$amp_acc <- 0

m00 <- update(mod0, .~., REML = TRUE)
pp <- exp(predict(m00, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0



### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m00, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 100"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe
##                     trt2 treatment_rewarded    IT_mm hive beeID amp_acc
## 1 Rewarded 1st 50 Buzzes                  F 4.267292    7 99999       0
## 2 Rewarded 2nd 50 Buzzes                  F 4.267292    7 99999       0
## 3 Rewarded 1st 50 Buzzes                  T 4.267292    7 99999       0
## 4 Rewarded 2nd 50 Buzzes                  T 4.267292    7 99999       0
##        blo      bhi predMean
## 1 49.00457 57.96004 53.73907
## 2 44.26685 52.01662 47.85855
## 3 42.39938 50.11380 46.54231
## 4 39.75753 48.64500 43.76588

Make amplitude plots

# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
facLevs = c("Rewarded 1st 50 Buzzes\n (buzz #51-100)", 
               "Rewarded 2nd 50 Buzzes\n (buzz #1-50)", 
               "Rewarded 1st 50 Buzzes\n (buzz #1-50)",
               "Rewarded 2nd 50 Buzzes\n (buzz #50-100)"
               )

pframe$t2 <- factor(facLevs,  levels = facLevs[c(3,1,2,4)])


g00 <- ggplot(pframe, aes(x=treatment_rewarded, color = treatment_rewarded, y=predMean))+
     geom_point()+ 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Bee Rewarded") + 
    theme(legend.position = "none", 
          legend.background = element_rect(fill=alpha('gray90', 0.5))) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded?", begin =0.3, end = 0.8)
g00

ggsave(plot = g00, filename = file.path(figDir, "StartStopPreds_amp2.pdf"), width =5, height = 3.5)


fac2Levs = c("51-100","1-50", "1-50","51-100")
              

pframe$buzzNums <- factor(fac2Levs,  levels = fac2Levs[c(3,1)])
pframe$beesRewarded = mapvalues(pframe$treatment_rewarded, from  = c(" F", " T"), 
                                to = c("No", "Yes"))

g001 <- ggplot(pframe, aes(x=buzzNums, color = beesRewarded, y=predMean))+
     geom_point()+ 
     labs(y =expression ("Sonication amplitude "(m~s^{-2})), x = "Sonication Number") + 
    theme(legend.position = "none", 
          legend.background = element_rect(fill=alpha('gray90', 0.5))) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded", begin =0.3, end = 0.8) +
  facet_wrap(~trt2)
g001

ggsave(plot = g001, filename = file.path(figDir, "StartStopPreds_amp3_mod0.pdf"), width =5, height = 3.5)


g001 <- ggplot(pframe, aes(x=beesRewarded, color = beesRewarded, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Sonication Number") + 
    theme(legend.position = c(.74,.8), 
          legend.background = element_rect(fill=alpha('gray90', 0.5))) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded", begin =0.3, end = 0.8) 
g001

bootstrap for a different model

# set number of bootstrap replicates for models
nbootSims = 5000

table(bees$hive) # more trials from hive 3
## 
##    7    8 
## 3887 1964
# using hive 7, since it's the one with the most data

pframe <- data.frame(trt2 = rep(levels(droplevels(bees$trt2)), 2),
                     treatment_rewarded = rep(levels(bees$treatment_rewarded), each = 2),
                     IT_mm = ITmean, 
                     hive = factor(7, levels = levels(bees$hive)),  
                     beeID = 99999)
pframe$amp_acc <- 0

m00 <- update(m2, .~., REML = TRUE)
pp <- exp(predict(m00, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0



### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m00, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 5000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe
##                     trt2 treatment_rewarded    IT_mm hive beeID amp_acc
## 1 Rewarded 1st 50 Buzzes                  F 4.267292    7 99999       0
## 2 Rewarded 2nd 50 Buzzes                  F 4.267292    7 99999       0
## 3 Rewarded 1st 50 Buzzes                  T 4.267292    7 99999       0
## 4 Rewarded 2nd 50 Buzzes                  T 4.267292    7 99999       0
##        blo      bhi predMean
## 1 47.85513 54.45403 51.03044
## 2 47.85513 54.45403 51.03044
## 3 42.39319 48.27959 45.20699
## 4 42.39319 48.27959 45.20699
pframe$beesRewarded = mapvalues(pframe$treatment_rewarded, from  = c(" F", " T"), 
                                to = c("No", "Yes"))

g001 <- ggplot(pframe, aes(x=beesRewarded, color = beesRewarded, y=predMean))+
     geom_point()+ 
     labs(y =expression ("Sonication amplitude "(m~s^{-2})), x = "Bee rewarded") + 
    theme(legend.position = "none", 
          legend.background = element_rect(fill=alpha('gray90', 0.5))) +
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) + 
  scale_color_viridis(discrete = TRUE, name  = "Bees Rewarded", begin =0.3, end = 0.8)
g001

ggsave(plot = g001, filename = file.path(figDir, "StartStopPreds_amp3_finalMod.pdf"), width =5, height = 3.5)

Amplitude GAMM4

GAMM4 shows time series

# start with gamm so I can show change by visit number
g01 = gamm4(log(amp_acc) ~ s(visitNum, by = trt2) + IT_mm + hive , random =  ~(1|beeID), data = bees)
par(mfrow = c(2,2))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, residuals = TRUE,pch = 20)

summary(g01$gam) # Summary for paper 
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(amp_acc) ~ s(visitNum, by = trt2) + IT_mm + hive
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.40245    0.26419   9.094  < 2e-16 ***
## IT_mm        0.34419    0.06234   5.521 3.51e-08 ***
## hive8       -0.23957    0.04982  -4.808 1.56e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                          edf Ref.df      F  p-value    
## s(visitNum):trt2Rewarded 1st 50 Buzzes 2.821  2.821 17.957 2.05e-10 ***
## s(visitNum):trt2Rewarded 2nd 50 Buzzes 2.150  2.150  3.144   0.0542 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0657   
## lmer.REML =  11311  Scale est. = 0.39028   n = 5851
summary(g01$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 11310.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5451 -0.4864  0.1642  0.6913  2.4609 
## 
## Random effects:
##  Groups   Name                                   Variance Std.Dev.
##  beeID    (Intercept)                            0.04546  0.2132  
##  Xr.0     s(visitNum):trt2Rewarded 2nd 50 Buzzes 0.01600  0.1265  
##  Xr       s(visitNum):trt2Rewarded 1st 50 Buzzes 0.02549  0.1596  
##  Residual                                        0.39028  0.6247  
## Number of obs: 5851, groups:  beeID, 96; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                2.402447   0.264191   9.094
## XIT_mm                                      0.344192   0.062339   5.521
## Xhive8                                     -0.239566   0.049822  -4.808
## Xs(visitNum):trt2Rewarded 1st 50 BuzzesFx1  0.026567   0.050437   0.527
## Xs(visitNum):trt2Rewarded 2nd 50 BuzzesFx1 -0.002741   0.046485  -0.059
## 
## Correlation of Fixed Effects:
##             X(Int) XIT_mm Xhive8 X(N15B
## XIT_mm      -0.993                     
## Xhive8       0.106 -0.183              
## X(N):2R150B  0.010 -0.010  0.012       
## X(N):2R250B  0.000  0.002  0.028  0.001
dev.off()
## null device 
##           1
# save gamm plots
pdf(file.path(figDir, "Gamm_amp_startStop_1st50.pdf"), width = 4.5, height = 3.5)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " log amplitude ", (m~s^{-2})))
plot.gam(select = 1, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded first", side = 2, padj = -3.5)
abline(v=50, lty = 2)
dev.off()
## null device 
##           1
pdf(file.path(figDir, "Gamm_amp_startStop_2nd50.pdf"), width = 4.5, height = 3.5)
par(mai= c(0.9,1,0.3,0.3))
dd = expression(paste('Estimated ', Delta, " log amplitude ", (m~s^{-2})))
plot.gam(select = 2, xlab = "Visit number", ylab = dd, g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE, bty = "l")
mtext("bees rewarded second", side = 2, padj = -3.5)
dev.off()
## null device 
##           1